Optimal Fourier filtering of a function that is strictly confined within a sphere 
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We present an alternative method to filter a distribution, that is strictly confined within a sphere 
of given radius r^, so that its Fourier transform is optimally confined within another sphere of 
radius kc- In electronic structure methods, it can be used to generate optimized pseudopotentials, 
pseudocore charge distributions, and pseudo atomic orbital basis sets. 

PACS numbers: 71.15.-m 



In some computational problems we are interested in 
distributions that are strictly confined within a sphere of 
given radius (i. e. defined to be strictly zero outside that 
sphere) and, simultaneously, optimally confined within 
another sphere in reciprocal space, so that they can be 
well approximated by a finite number of Fourier compo- 
nents or, equivalently, by a finite number of grid points 
in real space. Within the field of electronic structure 
calculations, this typically occurs in the real-space ap- 
plication of pseudopotentials ^"^ . In the specific case of 
the SIESTA density functional method^'^, this problem 
arises in the evaluation, using a real-space grid, of matrix 
elements involving strictly localized basis orbitals* and 
neutral-atom pseudopotentials. Those integrals produce 
an artificial rippling of the total energy, as a function 
of the atomic positions relative to the grid points (the 
so-called "egg box" effect), which complicates consider- 
ably the relaxation of the geometry and the evaluation 
of phonon frequencies by finite differences. 

We have proposed recently a method to filter a distri- 
bution simultaneously in real and reciprocal space^. Such 
filter is optimal, in the sense of minimizing the norm of 
the function outside two spheres of radius Tc and kc in 
real and Fourier space, respectively. It works by project- 
ing the distribution to be filtered on a basis of functions 
that have the same shape in real and reciprocal space, 
and that are thus optimally confined in both. However, 
because of the uncertainty principle, such basis functions, 
and the resulting filtered distribution, cannot be strictly 
confined in any of the two spaces. Thus, if we insist in the 
strict confinement in real space, and therefore we trun- 
cate the filtered pseudoatomic orbitals beyond r^ they 
will have a discontinuity at re, and therefore an infinite 
kinetic energy. In practice, the smallness of the discon- 
tinuity, and the use of integration grids with finite spac- 
ings, makes the problem more academic than real. But 
occasionally, when trying to converge the results to very 
high precision, it is annoying to have such a potential 
problem. Another, independent problem in our previ- 
ous procedure is that the resulting functions, that were 
expanded in Legendre polynomials, do not obey exactly 
the correct behavior for r — > 0. In the present work, we 
propose an alternative method to filter a distribution so 
that it is always strictly confined in real space, while it 
is optimally confined in reciprocal space. 



Consider an initial function with a well defined angular 
momentum and strictly confined within a sphere: 

p.(._(Fo{r)Yr{r) iir<rc 
°^ ' ~ \0 otherwise, ^ ' 

where Fo{r) is continuous and Fo{rc) = 0. We are using 
the same symbol for Fo{r) and its radial part FQ{r), since 
it does not lead to any confusion. Y™{r) is a real spher- 
ical harmonic. To require that F{r) [the filtered version 
of -Fo(r)] remains strictly zero for r > r^ and continuous 
at Tc (so that its kinetic energy is finite), we will expand 
it in terms of spherical Bessel functions ji with a zero at 
re- 

F(r) = I Sn=i '^nNinjiihnr) if r < r^ ^2) 

[0 otherwise, 

where M is large enough to represent the function with 
the required accuracy, kinTc is the nth root oi ji(x), and 
Nin are normalization constants given by 



N^n^ 



r'drjt{k,r,r)^^jt+i{kinrc). (3) 



The Fourier transform of F{v) is 



G(k) 



(27r)3/2 



d^r e-*'^""F(r) = G(fc)Y,"(k), (4) 



where we have introduced the factor i' to make G'(k) real, 
and 



N 



G(fc) = ^G„jta(A:), 



(5) 



n=l 



where jin{k) is the Fourier transform of Ninji{kinr): 
1 



jln{k) = 



(27r)3/- ,0 



47rr dr Ninji{hnr) ji{kr) 



cV^^ \ jl+liklnTc) 



if fc = kiji 



"" ^ ' k^ti^r- jiikrc) otherwisl^) 



The basis functions Ninji{kinr) are the solutions to 
Schrodinger's equation in a potential V(r) = for r < Vc- 



According to the variational principle, they are the func- 
tions that minimize the kinetic energy, among those 
strictly confined within a sphere of radius re- Some of 
the Fourier transforms jin{k) are shown in Fig. 1. They 
are delta-like functions in reciprocal space, broadened be- 
cause of their confinement in real space, according to the 
uncertainty principle. 



leads to the eigenvalue equation 




k/k 




k/k 



FIG. 1: Upper panel: squared Fourier transform, k^jin{k), 
of some normalized spherical Bessel functions ji (kinr) strictly 
confined in r < re. Lower panel: squared Fourier transform 
of selected solutions, gin{k), to the problem of minimizing the 
kinetic energy in fc > fcc, Eq. (7). I = 0, kcVc = 25. 



A conventional and straightforward method to filter 
F(r) would be to project it on the basis ji(fc(„r), with 
kin < kc, i. e. by truncating the series in Eq. (2). A better 
procedure is to use a basis of orthonormal functions that 
minimize not the total kinetic energy, but specifically the 
kinetic energy in the region k > kc that we want to filter 
out: 



k^dk gf{k) 



(7) 



Expanding the solutions in the primitive basis, 

9i{k)^^c.Jin{k), (8) 



A^ 



H„ 



ecr, 



(9) 



where e is a Lagrange multiplier to ensure normalization 
and 



Hn 



k'^dk jin{k) jim{k) 



ka 



k^ ^ 

"'ln"nrn 



k^dk 3in{k) ji„,{k). (10) 



The resulting eigenfunctions gin{k) (that we will call "fil- 
terets") are qualitatively very similar in real space to 
those in ref. [9] and therefore they are not reproduced here 
again. Fig. 1 shows them in reciprocal space for a very 
small value kcrc = 25, used to emphasize the effects of 
an extreme confinement. When fc;„ << fcc or fc(„ >> kc, 
they are similar to the primitive functions jin{k). For 
kin ^ kc, however, they are considerably better confined 
within k < kc- 

The eigenvalues e;„ of Eq. (9) give the integral of the 
kinetic energy "leaked" outside kc- 



k'dk gUk). 



(11) 



fcc 



As expected, these eigenvalues are e/„ ~ for fc;„ < kc 
and tin — kf^ for kin > kc- They are compared in Fig. 2 
with the same integral of the original functions and it can 
be seen that they are much smaller for kin < kc- 




FIG. 2: Leaked kinetic energy, above the cutoff kc, of the 
solutions gin{k) to Eq. (7) (eigenvalues ei„, Eq. (11), filled 
symbols) compared to the same integral f^ k'^dk ji„{k) for 
the confined spherical Bessel functions (empty symbols). I — 
0, kcrc = 25. 

Since the functions jf^ik) minimize the total kinetic 
energy, the decrease of kinetic energy in fc > fcc by gf^ik) 
must be at the expense of a larger increase in fc < kc, 
resulting in a net increase. To control this increase, we 
have found convenient to give a small weight (say w ~ 



0.1) to the kinetic energy in k < kc- This can be done 
simply by multiplying the last integral in Eq. 10 by a 
factor [1 — w). A very small value w = 10^^ was used 
in Fig. 1, just to break the degeneracy of the functions 
gin{k) with kin « kc- Larger values yield functions 
somewhat intermediate between both panels. 

The filtered function F(r) is then obtained by project- 
ing the original function Fo(r) over the subspace spanned 
by the "filterets" gf^ik) with a sufficiently low eigen- 
value (say ein/kf^ < 0.01). We have checked that the re- 
sulting scheme produces pseudoatomic orbitals, neutral- 
atom potentials'', and partial-core-correction densities^" 
that are free of the mentioned pathologies of the previous 
scheme^, and that reduce the "egg box" effect in SIESTA 
at least as well. Overall, however, the convergence tests 
yield rather similar results and therefore we do not repeat 
here the figures and tables of reference 9. 

Finally, a practical remark on the filtering procedure 
is appropriate. In grid-based methods^^, in which the ki- 
netic energy is calculated by finite differences, it is appro- 
priate to use a filtering cutoff kc given by the maximum 
plane wave vector that can be represented in the grid 
without aliasing"'^. In SIESTA, however, the dominant 
kinetic energy is calculated by well converged two-center 
integrals'' that do not contribute to the egg box effect. 
In this case, it is more convenient to fix kc by some in- 
dependent criterion, so that the orbitals (and the kinetic 
energy) do not depend on the integration grid used to 
calculate the exchange-correlation and pseudopotential 



interactions. Thus, we can fix an energy threshold ec 
such that 



ec = / k'^dk (j)^{k), 
'fee 



(12) 



where (j) ^re the original (unfiltered) atomic basis or- 
bitals. This criterion will yield different (but balanced) 
reciprocal-space cutoffs kc for each orbital, in the same 
spirit that the "energy shift" '' fixes their cutoffs Tc in 
real space. The grid cutoff will then be fixed to ^ 1.5 — 2 
times the maximum filter cutoff of all the orbitals (this 
factor coming from the fact that the plane wave cutoff for 
the density is larger than that for the wave functions). 

In conclusion, we have presented a simple but power- 
ful method to generate a basis of orthonormal functions 
("filterets"), with a given angular momentum, which are 
strictly confined within a cutoff radius in real space and 
optimally confined within another cutoff in Fourier space. 
We have described their use to filter a function that is 
strictly confined within a sphere. In addition, these or- 
thonormal functions constitute themselves a general and 
systematically improvable basis for converged calcula- 
tions using localized basis orbitals^^. This possibility will 
be explored in future works. 
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